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Abstract 

We investigate the excitonic energy transfer (EET) in the Fenna-Matthews-Olsen complex 
and obtain the linear absoiption spectrum (at 300 K) by a phenomenological time-convolutionless 
(TCL) master equation which is validated by utilizing Path Integral Monte Carlo (PIMC) sim- 
ulations. By applying Marcus' theory for choosing the proper Lindblad operators for the long- 
time incoherent hopping process and using local non-Markovian dephasing rates, our model 
shows very good agreement with the PIMC results for EET. It also correctly reproduces the 
linear absorption spectrum that is found in experiment, without using any fitting parameters. 

Introduction 

Since the seminal experiments on the Fenna-Matthews-Olsen (FMO) complex performed by the 
Fleming group in 2007- and subsequently confirmed by the Engel group in 2010, 2 many theoreti- 
cal models have been proposed in order to properly describe the observed long-lasting oscillations 
in the 2D spectroscopic data. While many methods are based on quantum master equations of 
some sort or the other, like Lindblad or Redfield equations,-"- only recently a more detailed de- 
scription of dissipative effects has been attempted in the form of Path Integral Monte Carlo (PIMC) 
calculations based on results from atomistic modeling combining molecular-dynamics (MD) sim- 
ulations with electronic structure calculations.- In contrast to a phenomenological modeling of the 
"environment" (protein scaffold, water molecules, etc.), the latter approach utilizes BChl-resolved 
spectral densities which can be directly incorporated in the PIMC calculations to obtain an exact 
account of the full exciton dynamics. 

Most quantum master equation approaches use special analytical forms of spectral densities, 
such as Ohmic or Lorentzian forms. Although this allows to obtain solutions for the equations, it 
remains to be a phenomenological ansatz. So far the results that have been obtained by hierarchical 
time-nonlocal master equations with a Lorentzian spectral density show reasonable agreement with 
the experimental results of both the absorption and 2D spectra at 77 K.— 

In contrast, numerically exact methods like PIMC simulations— or the QUAPI method—""— 
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have proven to be capable to produce the exact quantum dynamics of excitonic energy transfer 
over the experimentally relevant timescales. However, this comes at the price of rather large com- 
putational costs. Here, we combine the respective strengths of numerically exact and approxima- 
tive methods while overcoming their respective weaknesses: we use PIMC results for the exact 
quantum dynamics of the excitonic population dynamics over intermediate timescales (i.e. 600fs) 
which still allow for a fast production of the respective results, yet are sufficient to allow for a com- 
parison to a time-convolutionless (TCL) master equation based on a Lindblad approach; we further 
corroborate our results by comparison with PIMC data for up to 1 .5 ps. For a model dimer systems 
the authors have already successfully demonstrated such a concept. 16 Here, our ansatz captures the 
long-time behavior by relating to Marcus' theory of electron transport— ^2 because we assume that 
the long-time behavior is governed by a classical hopping process between the individual sites,— ^ 
see below. For the short-time dephasing behavior we introduce non-Markovian dephasing rates.— 
Since, in principle, we then obtain results for arbitrarily long times, we have an efficient yet very 
accurate way to calculate arbitrary transfer properties. Furthermore, we use the master equation 
to calculate the absorption spectrum of the FMO complex (at 300 K), an observable which can 
straightforwardly be obtained experimentally. 21 This opens the possibility to estimate the valid- 
ity of the underlying microscopic Hamiltonian as well as its respective parametrization based on 
recent results from mixed quantum-classical simulations 22 by getting into direct contact with ex- 
perimental data. 



Energy transfer on the FMO complex 

Microscopic description The dynamics of single excitations on the FMO complex is often de- 
scribed by a tight-binding Hamiltonian with 7 localized sites, corresponding to the 7 bacteri- 
ochlorophylls (BChls) of the FMO monomer.- 2 ^ 2 ^ The influence of the protein scaffold and solvent 
on the excitonic dynamics is treated, in the spirit of the Caldeira-Leggett model,— as a collection 
of harmonic modes that are linearly coupled to each BChl. Previous studies showed no significant 
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correlations between the bath induced energy fluctuations at different sites so we assume that 
each BChl is coupled to its own individual environment. The full Hamiltonian of the system can 
now be written as: 

H = H S + H B + H SB , (1) 

with 



Hs = ££„|n)(n| + H J mn \m){n\ , (2) 

H B = E + ^^A) > (3) 

Hsb = £|n)(n|(c ;! ^ + A,i cl) ) • (4) 

n 

The state \n) corresponds to the single-excitation state of site n, the parameter £„ denotes the energy 
gap between ground and excited state of site n, and describes the excitonic coupling between 
sites m and n. Furthermore, X nK , P nK , m nK and (O nK denote the position, momentum, mass and 
frequency of the bath oscillators, respectively. In the interaction Hamiltonian Hsb, the constants 
c nK (in units of eV/m) denote the coupling strength between site n and the bath modes. We have 

(cl) 

included the classical reorganization energies A„ ' as a counter- term in Hsb to prevent further 
renormalization of the site energies by the environment.-22i2£i2£ This quantity is defined as: 

A w = *r d(0 m, (5) 

where J n (co) (in units of 1/s) is the spectral density of the bath that is coupled to site n. In terms of 
the system parameters, it is given by: 

J >^ ) = ^2^^-^)- (6) 
n K zm nK w nK 

The precise numerical values of the different parameters entering in the expressions above were 
obtained from combined quantum-classical simulations for the full FMO complex including the 



solvent. 1 

Effective master equation approach 

We use now use the microscopic description of the FMO complex to set up a phenomenological 
second order time-local quantum master equation. In doing so, we are able to reproduce the dy- 
namics obtained from the PIMC simulations as well as extending it to, in principle, arbitrary long 
times. Additionally, our approach also allows to obtain results for the linear absorption spectrum 
which are in close accordance to experimental findings. 

The spectral density of the FMO complex— leads to reorganization energies A,^ of the order 
of 0.02 — 0.09 eV, which is comparable to the differences in the site energies £ n , while the exci- 
tonic couplings J mn are of the order of 1 meV. This implies that we can expect that the protein 
environment is relatively strongly coupled to the FMO complex and that it therefore leads to a 
strong damping for the population dynamics. This is also reflected by the results of the PIMC 
simulations. 9 

We now assume that in the long-time limit, after most of the coherences (i.e. off-diagonal el- 
ements of the reduced density matrix in the site-basis representation) in the system have decayed, 
EET can be described by a classical hopping process between the different sites (BChl's), that is 
induced by the protein environment. The transfer rates k mn have to satisfy detailed balance, ensur- 
ing a correct equilibrium state, and are assumed to follow from Fermi's golden rule. Furthermore, 

(c\) (c\) 

the rates should also depend on the reorganization energies A„ and A,„ of the baths that are 
coupled to the sites n and m, reflecting the differences in the coupling strengths of the protein en- 
vironment to each BChl. Unlike Forster theory, which assumes incoherent hopping between the 
energy eigenstates of H s ,— Marcus's theory of electron transport satisfies all these properties,— 1 ^ 
leading to transfer rates k mn of the form: 




(7) 
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with j8 = l/k B T and a£5 = a£ 1} +a! c1) . 

Aside from incoherent transfer between the sites, the environment also induces a strong dephas- 
ing on each site. In the framework of the second order TCL master equation, 20 these dephasing 
rates (in units of 1/fs) are given by: 

pt poo 

X n (t)=2Re ds \ do)J n (co) [coth (ft hco/ 2) cos(cqs) -isin((Os)] . (8) 
Jo Jo 

Here, we use the spectral densities J n (co) which have been obtained by MD simulations in Ref.— 
and numerically calculate the correlation function. 

The TCL master equation that describes the excitation dynamics can now be written as:— 

^ ee J2f(/)p(0 = ~[H s ,p{t)] + ®{t)p{t). (9) 

Our numerical results (not displayed) show that the Lamb shift term that usually appears in this 
equation, only leads to a negligible difference in both the population dynamics and the linear 
absorption spectrum (the position of the peak is shifted by approximately -1 meV). The dissipator 
@(t) is assumed to take the following Lindblad form, according to the considerations above: 20 

^(0P(0 = £ 7m ' lW ( L " ; ^ L »™ " \ { L lm L mn,pyj • (10) 

The Lindblad operators are defined by L mn = \m)(n\ and the rates by y mm (t) = X m {t) and y mn (t) = 
k mn for m 7^ n. The operators L mm model the dephasing process, while the operators L mn model 
the incoherent transfer between sites m and n. This choice of Lindblad operators will lead - in the 
long-time limit - to incoherent hopping transfer between the sites, which is different from, e.g., 
Redfield theory, which requires incoherent transfer between the eigenstates | yr) of Hs, leading to 
Lindblad operators of the form L ~ \ y n ){Wm\ 

However, we note that the equilibrium state of our master equation (p eq ) is slightly differ- 
ent from the one that follows from Marcus theory p eq ,dbi which is given by detailed balance, 
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lim,^oop„„(0 = (l/Z)exp(-j8e n ) and lim,-^ p mn (t) = for m ^ n, where p mn (t) = (m\p(t)\n). 
This can be shown by noting that @(t)p eq ,db = but [Hs,p e q,db] 0- This implies that p eq ^b is 
not a stationary state of our master equation. For the present calculation, the devations are only of 
the order of 1%, so we still expect our approach to give good results. 

Path Integral Monte Carlo simulations 

PIMC simulations allow to extract the exact quantum dynamics in the presence of a dissipative en- 
vironment, both for charge transport-^^-^ as well as energy transfer. 9 In short, the time evolution 
of the reduced density operator of a dissipative quantum systems is calculated by employing the 
path integral representation for the propagator according to the Feynman- Vernon theory— for fac- 
torizing or its extension to correlated initial preparations.— These path integrals are then evaluated 
by a stochastic sampling process based on Markov walks through the configuration space of all 
conceivable quantum paths which self-consistently emphasize the physically most relevant ones. 
While there is no limitation with respect to the choice of system parameters for which PIMC simu- 
lation are capable of producing numerically exact results, this approach is subject to the notorious 
'dynamical sign problem',— which reflects quantum-mechanical interferences between different 
system paths and results in an increase of the computational effort necessary to obtain statistically 
converged results which scales exponentially with the timescale over which the dynamics of the 
system is investigated. However, the presence of a dissipative environment substantially weakens 
these interference effects and therefore the sign problem. Furthermore, it allows for various effi- 
cient optimization schemes which lead to a further soothing of this computational bottleneck, thus 
significantly enlarging the accessible timescales.-^^y^ 

For the present case, we utilize the PIMC data presented in Ref . 9 to demonstrate the reliability 
of the master equation results and extend some of the former to longer timescales. To that extend, a 
factorizing initial preparation has been employed, where, resembling the situation prior to the cre- 
ation of an exciton, the bath modes initially are in thermal equilibrium with respect to themselves, 
while the exciton has been modeled to be either initially localized on one of the seven BChl sites 
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or in one of the seven excitonic eigenstates. 
Population dynamics 

In Fig. Q] we show the population dynamics that is obtained by solving the TCL master equation, 
Eq. ©, for initial conditions corresponding to a localization on the sites \n), e.g. p (0) = \n) (n\, in 
comparison to the corresponding numerically exact PIMC results. The dotted curves represents the 
latter and the solid lines represent the results from our master equation approach. Fig. |2]shows an 
extension of the results up to 1200 and 1500 fs for initial preparations in sites 1 and 6, respectively. 

In general we observe good quantitative agreement of our approach with the PIMC results 
for all localized initial preparations and over all observed timescales. The largest deviations are 
observed for an initial condition localized on site 4 for which the bath has the lowest reorganization 
energy (0.025 eV). Since our assumption of a classical hopping process at long-times is based on 
having a strong coupling to the environment, we would expect that our approach becomes worse 
with decreasing reorganization energy. Also, from Fig. |2]one observes a good agreement in the 
approach to equilibrium, although the decay is slightly slower than predicted by the PIMC results. 

Figure |3] corroborates our results. Here, the excitonic excitation is initially in one of the seven 
eigenstates \\ff n ) of Hs- Again we find very good agreement with the PIMC results, where once 
more the strongest deviations occur for the initial preparation exhibiting the largest population 
on site 4. We note that there is no fitting parameter involved. Introducing a parameter which 
interpolates between purely coherent and purely incoherent transfer, as in,— could lead to a 
further improvement of the agreement. Nevertheless, already this rather simple phenomenological 
model captures most of the details which are present in the PIMC calculations. Additionally, it 
allows for a computationally cheap calculation of the linear absorption spectrum. 
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Figure 1 : Comparison of population dynamics of the 7 different sites of the FMO complex obtained 
from the numerically exact PIMC results (circles) with the results from the TCL master equation 
approach (solid lines) for different localized initial conditions on the sites \n), n = 1 , . . . , 7. Note 
that the statistical error of the PIMC calculations is typically smaller than the symbol size. There- 
fore we do not show the error bars. 
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Table 1: The numerical values for the x-, y-, and z- component as well as the absolute value 
squared of the transition dipole moments jl m = (n m ,xiHm,y,Hm,z) m units of Debye [D], The z 
axis is chosen along the C3 -symmetry axis of the FMO complex, and the y axis is chosen to be 
parallel to the N B - N D axis of BChl 1. 



m 



1 



5 


6 


7 


-6.39 


5.16 


0.0 


0.0 


2.98 


-1.14 


-0.45 


2.29 


5.85 


41.09 


40.70 


35.52 



0.0 -6.10 

1.86 1.08 

6.07 1.66 

40.32 41.09 



-5.27 
-3.04 
-2.10 
41.47 



0.0 
2.49 
5.85 
40.45 



10 




11 



Linear absorption spectrum 

Due to the strong dephasing that is induced by the protein environment, the coherences (i.e., the 
off-diagonal elements of the reduced density matrix in the site basis) die out rather fast. Hence, 
the dynamics of the populations is mostly insensitive to the exact behavior of the coherences. 
However, other observables, such as the linear absorption spectrum, are very sensitive to the short- 
time behavior, where the coherences are still present. Being able to access the full dynamics of the 
reduced density matrix for, in principle, arbitrary timescale, now allows us to access in particular 
the linear absorption spectrum for the FMO complex. We calculate the linear absorption spectrum 
and compare it to the absorption spectrum which was measured in experiment— and to the one 
which has been computed with mixed quantum-classical simulations by Olbrich et ah— 

The linear absorption spectrum is given by the Fourier transform of the two-time correlation 
function of the transition dipole moment (TDM) operator fi\— 



where Jl{t) = e lHt ' h fl,e~ lHt ' h is the TDM operator in the Heisenberg picture and JX =£,„ \l m (|m)(0| + |0)(m|), 
with jl m the TDM vector of site m. The two-time correlation function is evaluated in the excitonic 
ground state W° = |0) (0| <8> ps,i-6-, with no excitations present. 

To compute this correlation function with our approach, we use the following expression:— 




(11) 




(12) 



where the vector-operator V(t) satisfies the TCL master equation 



dV(t) 



(13) 



dt 
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Figure 4: Comparison of the linear absorption spectra, computed with the TCL master equation 
(solid red), mixed quantum-classical calculations 22 (dashed green) and obtained by experiment— 
(dotted blue). All three spectra are overlaid such that the position of the peaks are shifted to that of 
the experimental result. The inset shows the spectra at their original positions. 
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with the initial condition 



V(0) = Ahr fi { W °}=n\0) (0| = 2>„ \m) (0| , (14) 

m 

where tr$ and tr^ denote the trace over the excitonic and environmental degrees of freedom, re- 
spectively. Due to the time dependence of the dephasing rates that enter into the generator Jz?(f), 
we cannot obtain an analytical solution for the linear absorption spectrum. Therefore, we solve 
Eq. (IT3T ) numerically and obtain the absorption spectrum with Eq. (PT21) . 

The numerical values of the TDM vectors p, m are obtained by combining the data from Refs.— 
and — together with their relative orientations with respect to the C3 symmetry axis, taken from 
Ref.— In tabled] we provide the computed values of the jU m 's. 

In Fig. © we show the numerical results for the linear absorption spectrum and compare it 
to the experimental results by Freiberg et. ah— and computational results obtained by Olbrich 
et. at— utilizing the same parametrization of the Hamiltonian, Eq. CO, but based on ensemble- 
averaged wave packet dynamics.— We observe that both the spectra from the latter and from 
our theoretical approach are shifted by —0.04 eV and —0.02 eV, respectively, compared to the 
experimental spectrum. When overlaying those spectra such that the position of the peaks match, 
we observe almost perfect agreement for the overall actual line shape, even though our method 
lacks the detailed features of the experiment. We note that our approach yields a shift from the 
experimental findings which is only half as large as the one presented in Ref.— This indicates 
that despite the strong dissipation, quantum effects still play an important role for the excitonic 
dynamics which are not captured by the ensemble-averaged wave packet dynamics. 

Summary 

We have presented an approach to calculate the EET as well as the absorption spectrum of the 
FMO complex. Our approach is based on a phenomenological TCL master equation: The dis- 
sipator in the master equation is determined by incoherent hopping rates obtained from Marcus' 
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theory as well as by non-Markovian (pure) dephasing rates obtained from the bath autocorrelation 
function. Here, we have used spectral densities calculated from atomistic simulations of Ref.— We 
have demonstrated the quantitative reliability of our TCL master equation by a comparison with 
population dynamics obtained from numerically exact PEVIC simulations exhibiting an excellent 
accuracy for timescales up to the picosecond range, both, for the exciton initially found localized 
on any BChl as well as in an excitonic eigenstate. Furthermore, we also found very good agree- 
ment with the experimentally measured absorption spectrum. Here, the overall line shape of the 
absorption spectrum is reproduced very nicely, although the position of the peak is shifted by 0.02 
eV as compared to the experimental result, yet considerably less than for calculations based on 
wave packet dynamics. Since we have correctly reproduced both the population dynamics and the 
dynamics of the coherences, as reflected in the absorption spectrum, we also expect to be able to 
reproduce the results for 2D electronic spectroscopy that were found in experiments at ambient 
conditions, which is subject of current research. 
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